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SEMIANNUAL  TECHNICAL  REPORT  SUMMARY 


Contract  Objectives 

A  major  consideration  in  the  design  of  excavations  in  rock  is 
the  evaluation  of  the  structural  stability  of  the  opening.  An 
essential  step  in  this  evaluation  is  the  determination  of  the 
mechanical  state  in  the  rock  mass  in  the  vicinity  of  an  opening 
through  the  formulation  and  solution  of  boundary  value  problems. 
The  finite  element  method  has  proven  to  be  a  very  powerful  tool 
in  the  solution  of  boundary  value  problems.  Within  the  general 
framework  of  the  finite  element  method,  procedures  have  been 
developed  to  model  the  response  of  certain  classes  of  joints, 
cracks  and  fissures  in  the  rock  mass,  the  inability  of  rocks  to 
withstand  tension,  localized  yielding  of  rock  due  to  stress  con¬ 
centrations  and  the  time-dependent  (creep)  response  of  rock. 
While  preliminary  analysis  has  indicated  that  these  procedures 
have  the  potential  for  predicting  performance,  their  use  in 
design  will  be  limited  unless  their  reliability  is  established. 
The  objective  of  this  contract  is  to  evaluate  the  ability  of 
certain  available  finite  element  techniques  for  the  solution  of 
plane  problems  to  predict  the  stresses,  strains  and  deformations 
in  a  rock  mass  surrounding  an  excavation. 

General  Approach  and  Technical  Results 

The  approach  to  this  study  can  be  divided  into  two  phases:  (1) 
Developing  computational  procedures,  which  incorporates  the 
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essential  features  known  to  be  present  in  a  rock  mass,  for  de¬ 
termining  its  mechanical  state,  and  (2)  Analysis  of  case  histories 
and  comparing  predicted  and  measured  performance.  This  report 
is  concerned  primarily  with  Phase  1. 

The  essential  features  of  the  various  finite  element  techniques 
to  conduct  the  following  analyses  were  reviewed:  (I)  "No-tension" 
analysis,  (II)  Joint  Perturbation  Analysis,  (III)  Elastic-Plastic 
Analysis,  and  (TV)  Time-dependent  (Creep)  Analysis.  Because  the 
time-dependent  (creep)  response  of  the  great  majority  of  rocks 
is  not  as  significant  as  the  time -independent  response,  and  fur¬ 
thermore,  since  the  computational  procedures  for  time -independent 
and  time -dependent  problems  are  significantly  different,  it  was 
decided  that  only  time-independent  analyses  for  plane  problems 
would  be  considered  at  the  present  time.  It  was  found  that  dif¬ 
ferent  computational  procedures  were  utilized  to  conduct  the 
above-mentioned  time-independent  analyses.  In  order  to  develop 
a  single  computer  program  which  could  be  used  for  conducting 
analyses  which  include  the  features  in  (I),  (II),  and  (III),  it 
was  necessary  to  formulate  a  consistent  computational  procedure 
and  develop  a  program  on  the  basis  of  such  a  procedure.  The 
"initial  stress"  (stress  transfer)  technique  presented  by  Zienkiewicz 
and  his  co-workers  was  used  and  an  operational  program  which  had 
the  capability  of  conducting  a  "no  tension,"  joint  perturbation 
and  elastic-plastic  analysis  was  developed.  Various  illustrative 
examples  were  solved  using  the  combined  program  and  the  results 
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compared  when  possible  with  results  published  by  other  investi¬ 
gators.  Except  for  the  case  of  an  elasto-plastic  analysis  for 
a  material  with  frictional  characteristics,  where  convergence 
difficulties  were  encountered,  the  results  obtained  appeared 
satisfactory.  Further  work  is  being  done  to  minimize  conver¬ 
gence  difficulties  in  the  elasto-plastic  analysis. 

POD  Implications 

The  evaluation  of  the  structural  stability  of  underground  openings 
is  an  essential  step  in  design  of  underground  structures  in  rock. 
The  use  of  such  structures  is  likely  to  increase  in  future 
for  civilian  and  military  purposes.  The  development  of  a  single 
finite  element  program  which  includes  the  capability  of  realis¬ 
tically  modelling  rock  characteristics  will  be  of  considerable 
assistance  in  the  design  of  underground  structures. 

Considerations  for  Further  Research 

The  next  phase  of  the  program  deals  with  establishing  an  estimate 
of  the  ability  of  the  computational  procedures  in  predicting 
stresses  and  deformation  in  rock  masses.  Further  research  should 
consider  increasing  the  capability  of  the  computational  proce¬ 
dures  to  include  construction  sequence  and  rock  reinforcement 
schemes. 
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INTRODUCTION 

The  development  of  theoretically  sound  methods  for  designing 
excavations  in  rock  is  of  significance  because  of  the  increased 
use  of  underground  facilities  for  civilian  and  military  purposes. 

A  major  consideration  in  the  design  of  excavations  in  rock  is 
the  evaluation  of  the  structural  stability  of  the  opening.  If 
a  mechanistic  approach  is  taken  for  the  evaluation  of  the  struc¬ 
tural  stability,  then  an  essential  step  in  the  evaluation  is 
the  determination  of  the  mechanical  state  (i.e.  stresses,  strains, 
and  deformation)  in  the  rock  mass  in  the  vicinity  of  the  excava¬ 
tion.  In  recent  years  new  technology  has  been  developed  which 
has  the  potential  of  predicting  with  greater  accuracy,  than 
heretofore  possible,  the  mechanical  state  in  the  rock  mass. 

This  would  place  the  evaluation  of  the  structural  stability 
of  an  opening  in  rock  on  a  theoretically  sound  basis.  This 
new  technology  has  been  primarily  in  the  area  of  numerical 
methods  for  the  solution  of  boundary  value  problems.  However, 
the  application  of  these  methods  to  practical  design  problems 
has  been  very  limited.  The  primary  reason  for  this  is  the 
designer's  lack  of  confidence  in  the  ability  of  theoretical 
methods  to  assist  him  in  analyzing  practical  design  problems. 

Of  the  numerical  methods  developed,  the  finite  element  method 
has  proven  to  be  the  most  powerful  for  the  solution  of  boundary 
value  problems,  (Clough,  1960,  Wilson,  1963,  1965,  Zienkiewicz, 
1967)  .  The  initial  application  of  the  finite  element  method  to 
rock  mechanics  problems  used  techniques  (computer  programs)  that 
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had  been  developed  for  the  structural  analysis  of  linear 
elastic  continuous  structures.  When  the  results  obtained 
from  using  these  techniques  were  compared  with  field  results, 
they  were  found  to  be  inadequate  for  predicting  the  response 
of  rock  masses  in  the  vicinity  of  an  excavation.  Typical 
results  from  such  studies  have  been  presented  by  Judd  and 
Perloff  (1971).  It  was  recognized  that  one  of  the  major 
reasons  for  the  discrepancies  between  computed  and  observed 
behavior  was  the  inability  of  the  techniques  utilized  to  make 
the  computations  to  include  the  natural  geologic  discontinuities 
that  were  usually  prevalent  in  a  rock  mass.  These  discontinuities 
could  exist  in  the  rock  formation  prior  to  the  excavation  or  could 
be  a  result  of  the  induced  stresses  due  to  the  excavation  causing 
failure  in  the  rock.  Methods  of  analysis  which  have  the  capa¬ 
bility  of  modelling  joints  and  other  forms  of  discontinuities 
that  commonly  exist  in  rock  masses  have  been  developed  by 
Goodman,  Taylor,  and  Brekke  (1968),  and  Zienkiewicz,  et  al 
(1970).  Methods  for  stress  analysis  in  a  rock  mass  which  can¬ 
not  sustain  tensions  due  to  the  presence  of  cracks  and  fissures 
have  also  been  developed,  Zienkiewicz,  et  al  (1968).  In 
addition,  methods  for  analyzing  localized  failure  in  a  rock 
mass  due  to  yielding  have  been  presented  by  Reyes  and  Deere 
(1966) .  Approximate  techniques  for  incorporating  nonlinear 
time-dependent  material  properties  have  been  developed  by  Nair 
and  Boresi  (1970) .  Preliminary  analysis  has  indicated  that 
these  new  methods  of  analysis  using  finite  element  techniques 
have  the  potential  of  predicting  performance  with  improved 
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accuracy.  However,  th^  methods  of  analysis  which  account  for  a 

1  ‘ 

specific  rock  characteristic  e.g.,  no  tension  or  joints,  are 

’  III 

not  adequate  to  study  the  general  case  where  all  these  factors  may 

,  i  '  (  i  '  i  i  1 

be  present.  Furthermore,  there  has  been  very  limited  verification 
of  thbse  techniques  on  the  basis  of  comparison  with  measured 

1  i  '  '  i  ‘ 

field  performance.  Without  field  verification,  the  use  of  these 

i  1  '  i 
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new  techniques  in  the  design  of  excavations  in  rock  will  remain 

limited,  Therefore,,  to  develop  a  theoretically  sound  method  1 
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for  designing  excavations  in  rock,1  which  will  be  used  in  prac- 
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tice,  an  essential  step  is  to,  establish  the  reliability  of  the 
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available  methods  of  stress  analysis  in  predicting  the  behavior  , 

i  '  ,  i 

of  rock  masses.  • 
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PURPOSE  " 
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The  purpose  of  thisi  study  is  to  evaluate  the  reliability  of  the 

1  ■  i  i  ,  .*  1 

available  finite  element  techniques  for  the  solution  of  plane  1 

i  1  i 

problems  ,to  predict  the  stresses,,  strains,  and  deformations  in 

,  i  , ,  1  i  .  i  i  1  i 

a  rock  mass  surrounding  an  excavation. 
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GENERAL  APPROACH  11 
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t  . 

The  approach  to  this  study  can  be  divided  into  two  phalses:  (1) 

1  i  ,  i  1  i 

Developing  a  general  analytical  procedure  (i.e.,1  a  finite  ele-  ■ 
ment  computer  program  with  Wide  capabilities)  for  determining 

"«  ,  '  1  i 

the' mechanical  state  in  a  rock  mass,  and  (2)  Analysis  of  case 
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histories  to  compare  predicted  and  observed  values,  i 

1  i  i 
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SCOPE  OF,  THIS  REPORT*  .  ,  .  , 

This  report  deals  primarily  with  Phase  1.  Certain  preliminary 

1,1!  1 

work  connected  with  Phase  2  is  summarized, 
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Phase  I  consisted  of  three  steps: 
i  a.  Review  of  pertinent  available  finite  element  programs. 

b.  Consolidation  of  the  various  programs  into  a  single 

I 

general  program. 

i  i 

I 

c.  Use  of  examples  to  illustrate  the  use  of  the  program. 


Preliminary  work  has  been  conducted  on  the  selection  of  case 

i  1 

histories  for  analysis  and  the  availability  of  three  dimensional 

I 

finite  element  programs.  ( 

I  i 
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ESSENTIAL  FEATURES  OF  FINITE  ELEMENT  TECHNIQUES  FOR  MODELLING 
ROCK  BEHAVIOR  (Phase  la) 

,  I 

Various  techniques  have  been  utilized  in  conjunction  with  the 
finite i element  method  to  include  nonlinear  and  time -dependent 

I 

properties,  elasto-plastic  yielding,  and  geologic  discontinuties 

1  i 

in  the  stress  analysis  of  rock  masses.  These  techniques  were 

( 

i  !  1  1  i 

evaluated  for  the  purpose  of  developing  a  general  computer  pro- 

i  i  ' 

gram  for  plane  problems  which  can  incorporate  the  significant 

i 

features  that  exi^t  in  rock  masses.  The  available  techniques 
use  different ■ computational  techniques  within  the  general  frame¬ 
work  of  the  finite  element  method  to  develop  solutions  to 

i  (  r  f 

boundary  value  problems.  The  basic  concepts  of  the  finite 
element  method  have  been  discussed  extensively  in  the  litera¬ 
ture,  e.g.  Clough  (1960) ,  Wilson  (1963,  1965)  and  Zienkiewicz 

I 

(1967).  These  will  not  be  repeated  here.  Those  aspects  which 

i  i 

! 

are  considered  special  modifications  for  rock  mechanics  problems 
are  discussed  in  general  terms  with  respect  to  the  computational 

I  !  i  I 

techniques'  in  the  subsequent  paragraphs.  These  are  considered 

•  i 

i  5 
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in  four  categories  (1)  "No  Tension"  analysis,  (2)  Joint  Perturba¬ 
tion  Analysis,  (3)  Elasto-Plastic  Analysis,  and  (4)  Time  dependent 
analysis.  Detailed  explanation  of  these  techniques  is  provided 
in  a  subsequent  section  where  the  basis  of  the  developed  computer 
program  is  described. 

1 .  "No  Tension"  Analysis 

When  numerous  cracks  and  fissures  are  present  in  a  rock  mass,  it  has 
been  assumed  that  the  rock  is  incapable  of  withstanding  tensile 
stresses.  A  procedure  for  modelling  this  nonlinear  behavior  has  been 
presented  by  Zienkiewicz,  et  al  (1968).  This  method  is  called  "no 
tension"  or  "stress  transfer"  analysis.  This  method  is  composed 
of  the  following  four  essential  steps: 

(a)  A  linear  elastic  solution  to  the  problem  is  first 
obtained.  The  induced  changes  in  stress  are  added 
to  the  initial  stresses  and  the  principal  stresses 
computed. 

(b)  Those  elements  in  which  tensile  stresses  are  present 
are  identified.  As  the  material  is  assumed  incapable 
of  sustaining  them,  the  calculated  tensile  principal 
stresses  are  eliminated  without  permitting  any  point 
in  the  structure  to  displace.  In  order  to  maintain 
equilibrium,  equivalent  (balancing)  nodal  forces 

are  calculated  and  temporarily  applied  to  the  structure. 

(c)  The  elastic  analysis  is  repeated  (i.e.,  the  equilibrium 
equations  are  satisfied)  to  remove  the  balancing  nodal 
forces  and  the  check  for  tensile  stresses  is  repeated. 
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(d)  If  at  the  end  of  step  (c)  tensile  stresses  are  still 
present,  steps  (b)  and  (c)  are  repeated  until  no 
appreciable  difference  in  magnitude  and  distribution 
of  stresses  is  obtained  upon  further  iteration.  It 
•  —is  important  to  recognize  that  the  initial  elastic 
stiffness  matrix  is  used  throughout  the  computation. 

This  type  of  analysis  has  been  used  by  Zienkiewicz,  et  al  (1968) 
in  the  solution  of  plane  strain  problems  in  rock  mechanics.  Heuze, 
Goodman,  and  Bornstein  (1971)  have  used  the  same  technique  to  study 
the  deformability  of  jointed  rock  in  borehole  jacking  tests.  The 
convergence  of  the  solution  has  been  found  to  be  very  slow.  The 
number  of  iterations  required  for  convergence  ranged  from  10  to 
15  cycles.  Zienkiewicz,  et  al  (1968)  reported  that  faster  con¬ 
vergence  could  be  obtained  if  the  elastic  constants  are  modified 
by  reducing  the  modulus  in  the  direction  of  the  principal  tensile 
stress.  With  this  modification  the  stiffness  matrix  would  have  to 
be  recomputed  at  every  stage. 

2 .  "Joint  Perturbation"  Analysis 

A  one-dimensional  joint  element  was  developed  by  Goodman,  Taylor, 
and  Brekke  (1968)  and  applied  to  several  rock  mechanics  problems. 

This  technique  is  capable  of  modelling  the  behavior  of  joints, 
bedding  planes  and  other  geologic  discontinuities.  The  original 
version  defined  normal  and  shear  stiffnesses,  and  Kg,  for  joint 
elements  and  incorporated  the  stiffness  of  these  elements  into 
the  stiffness  of  the  overall  structure.  Fig.  1  shows  a  one-dimensional 
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FIG.  1  -  LINKAGE  OR  "JOINT"  ELEMENT  WITH 
ITS  LOCAL  COORDINATE  SYSTEM 


joint  element  with  its  local  coordinate  system  and  the  sign 
convention  on  stresses  and  displacements.  The  joint  element 
was  assumed  to  have  no  thickness.  The  original  version  was 
only  valid  for  a  linear  elastic  analysis.  In  his  study  on 
the  deformability  of  joints,  Goodman  (1969)  has  reported  that 
there  are  four  typical  shear  stress-deformation  relationships 
for  various  weak  surfaces  as  shown  in  Fig.  2.  It  is  apparent 
that  peak  shear  strength  in  most  cases  is  greater  than  residual 
shear  strength.  After  the  peak  shear  strength  is  reached,  shear 
stress  drops  to  residual  values  and  appreciable  movement  can 
take  place  without  increase  in  shear  stress.  It  is  obvious 
that  this  type  of  stress -deformation  behavior  cannot  be  approxi¬ 
mated  by  one  single  value  of  shear  stiffness  as  used  in  the 
original  analysis,  Goodman  et  al  (1968).  Heuze,  Goodman  and 
Bornstein  (1971)  modified  the  original  version  so  that  normal 
and  tangential  properties  can  be  varied  as  a  function  of  dis¬ 
placement.  They  used  an  iterative  approach  to  solve  this  type 
of  nonlinear  problems.  The  same  boundary  value  problem  is 
analyzed  repeatedly.  Each  time  a  new  set  of  normal  and  tan¬ 
gential  stresses  and  disj lacements  across  a  joint  element  is 
calculated  and  used  together  with  joint  properties  and  joint 
strength  to  modify  the  normal  and  shear  stiffness  for  the  next 
iteration.  In  this  approach  a  new  stiffness  matrix  has  to  be 
formulated  and  computed  for  every  iteration.  This  is  a  signifi¬ 
cant  computational  difference  from  the  stress  transfer  technique 
used  in  the  no  tension  analysis  where  the  stiffness  matrix  is 
not  changed  with  repeated  iterations. 
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Zienkiewicz,  et  al  (1970)  have  shown  that  an  iterative  process 
similar  to  "stress  transfer"  analysis  may  be  employed  in  the 
analysis  of  jointed  rock  systems. 

3.  Elasto-plastic  Analysis 

An  elasto-plastic  analysis  has  been  suggested  by  various  re¬ 
searchers  to  account  for  the  possible  yielding  of  rock  due 
to  the  stress  concentrations . induced  in  the  rock  mass  around 
an  underground  opening.  Reyes  and  Deere  (1966)  developed  a 
method  based  on  the  incremental  theory  of  plasticity  to  study 
elasto-plastic  behavior  of  underground  openings.  From  the 
computational  point  of  view  the  process  used  by  Reyes  and 
Deere  (1966)  has  the  disadvantage  in  that,  at  each  iteration 
the  stiffness  of  the  structure  is  changed,  requiring  a  re¬ 
formulation  and  computation  of  the  stiffness  matrix  which 
will  involve  extra  computer  time.  An  alternative  approach 
has  been  developed  by  Zienkiewicz,  et  al  (1969).  This  approach 
uses  a  technique  referred  to  as  the  "initial  stress"  technique 
which  is  consistent  with  the  "no  tension"  analysis.  The  solu¬ 
tion  obtained  by  the  "initial  stress"  process  involves  a  series 
of  load  increments.  For  each  load  increment,  the  solution 
satisfying  equilibrium  and  yield  criteria  is  achieved  by  a 
series  of  approximations  or  iterations.  Fig.  3  illustrates 
how  the  final  solution  may  be  achieved  by  a  series  of  iterations. 
Baker,  Sandhu,  and  Shieh  (1969)  have  also  developed  a  computa¬ 
tional  technique  for  the  elasto-plastic  analysis.  The  basic 
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concept  is  quite  similar  to  that  used  in  the  initial  stress 
approach  presented  by  Zienkiewicz,  et  al  (1969).  To  save 
computer  execution  time  both  approaches  use  the  initial  (elas¬ 
tic)  stiffness  of  the  structure  at  each  step  of  computation. 
Corrective  body  forces  are  calculated  and  applied  to  the 
structure  during  iterations  to  insure  that  the  element  is 
just  on  the  yield  surface.  The  only  difference  between  the 
approaches  of  Zienkiewicz,  et  al  and  Baker,  et  al  is  that 
the  latter  employed  a  different  method  to  insure  the  convergence 
of  the  final  solution  for  each  load  increment. 

4.  Time -dependent  Analysis 

Time -dependent  analysis  for  problems  in  rock  mechanics  can 
be  considered  in  two  categories.  The  first  includes  those 
cases  where  the  boundary  value  problem  changes  with  time. 

Such  problems  include  consideration  of  the  gradual  (time- 
dependent)  creation  of  the  excavation  or  the  installation 
of  reinforcement  (e.g.  rock  bolts)  at  various  times  during 
the  construction  of  the  opening  or  later  during  its  performance. 

The  second  category  includes  those  problems  where  the  properties 
of  the  rock  surrounding  the  excavation  are  time-dependent  (creep). 
This  review  was  confined  to  problems  in  the  second  category,  and 
in  the  subsequent  discussions,  time-dependent  analyses  refer 
only  to  time -dependent  material  properties. 

Not  much  emphasis  has  been  placed  on  developing  time -dependent 
analyses  for  rock  mechanics  problems  because  the  great  majority 
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of  rocks  do  not  exhibit  significant  time -dependent  behavior. 

In  structural  analysis  two  approaches  have  been  taken  in  de¬ 
veloping  solutions  for  time -dependent  problems.  These  depend 
on  the  complexity  of  material  response.  For  linear  problems 
the  theory  of  linear  viscoelasticity  has  been  utilized.  For 
materials  exhibiting  nonlinear  behavior,  empirical  stress-strain 
relations  together  with  iterative  solution  techniques  have 
been  used.  Excessive  computer  costs  and  difficulties  asso¬ 
ciated  with  determining  reliable  material  properties  have 
limited  the  application  of  such  analytical  techniques  to 
problems  in  rock  mechanics. 

Time -dependent  aialyses  for  nonlinear  material  properties 
have  been  developed  by  Deere  and  Boresi  (1963),  Nair  (1967), 
Aiyer  (1969)  and  Nair  and  Boresi  (1970).  All  these  analyses 
were  developed  for  specialized  problems  and  with  the  exception 
of  Aiyer' s  work  were  confined  to  spherically  symmetric  or 
axisymmetric  problems.  The  approach  utilized  by  Nair  and 
Boresi  (1970)  uses  the  finite  element  method  and  uses  a  compu¬ 
tational  technique  which  can  be  applied  to  the  analysis  of 
plane  problems.  The  basic  concept  used  in  this  analysis  is 
similar  to  that  proposed  by  Greenbaum  (1966) .  The  method  of 
incorporating  the  time-dependent  behavior  of  the  material  is 
briefly  described  as  follows:  As  a  first  step  in  the  compu¬ 
tation,  an  elastic  stress  distribution  is  computed;  the  effec¬ 
tive  stress  is  then  computed.  Based  on  the  value  of  the  effec¬ 
tive  stress,  the  effective  strain  rate  is  computed  using  the 
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creep  stress-strain  time  law.  Taking  a  small  time  increment  At, 
the  incremental  strains  can  now  be  computed.  From  these  incre¬ 
mental  strains,  the  incremental  stresses  are  computed  using  the 
elastic  stress -strain  relations.  These  incremental  stresses 
are  then  input  into  the  program  as  initial  stresses,  converted 
into  nodal  point  forces,  and  the  problem  is  solved  as  an  elas¬ 
tic  problem  and  the  new  stresses  and  displacements  are  determined. 
These  stresses  are  then  used  to  develop  the  new  incremental 
strains  and  the  process  repeated.  In  this  manner,  the  variation 
of  stress  and  displacement  with  time  is  determined. 

Summary  of  Achievement  for  Phase  la 

A  review  of  the  available  literature  indicated  that  methods  are 
available  for  modelling  geologic  discontinuities  and  failure 
characteristics  that  occur  in  rock  masses.  However,  the  compu¬ 
tational  techniques  associated  with  the  various  methods  are 
not  identical.  Consequently,  the  first  step  in  the  development 
of  any  general  purpose  computer  program, which  would  include  the 
capabilities  of  the  individual  programs,  would  be  to  make  the 
various  computational  techniques  consistent. 

It  was  also  found  that  time -dependent  analyses  were  not  developed 
to  the  level  of  sophistication  of  the  time -independent  analyses. 
Furthermore,  the  cost  associated  with  a  time -dependent  analysis 
and  the  fact  that  relatively  few  rocks  exhibit  significant  time- 
dependent  behavior  leads  to  the  conclusion*  at  this  time,  that, 
the  development  of  a  general  purpose  program  for  the  time -independent 
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plahe  problems  be  treated  separately  from,  the  development  of  a 
program  for  time-dependent  analyses.  , 


.  i  i 

tJEVfiLOPMENT  OF  A  GENERAL  COMPUTER  PROGRAM  (Phase  lb)  '  , 

t  i  1 

In  order  to  develop  a  single  general  computer  program,  it'  was 

i  i 

i  necessary  to  develop  a  consistent  computational  technique  to  model 

(  1  /  if’ 

the  different  aspects  of  rock  behavior.  On  the  basis  of  the  ! 
review  of  the  available  techniques,  it  was  concluded  that  the 
"initial  stress"  (stress  transfer)  technique  presented  by  Zienkie- 

, '  1  ■  '  i  1  , 

wicz  and  his  co-workers  would  provide  a  consistent  approach  in 

;  I  f  '  !  ! 

the  development  of  a  consolidated  computer  program,  for  the 

1  1  i  i  1  f  i 

stress  analysis  of  plane  problems,  which  is  cdpable  of  modelling 

.  ■  t  •'  ' 

joints  and  other  geologic  discontinuities,  and  elasto-plastic 

■  ii  i 

or  time -dependent  rock  properties.  The  major  reason  for  selecting  1 

-  ,■  ’  I 

this  approach  was  the  computational  advantage  that  results  from 

1  '  I 

•  ■ 

using, the  initial  elastic  stiffness  at  every  stagd  of  the  solution 

1  1  1  *  f  1  '■  i  ,  , 

1  process.  The  features  incorporated  in  the  combined  comiputer  program 
for  time -independent  plane  problems  are:  ,  (1)  No  Tension  analysis, 

1  i  1 

i  1 

(2)  Joint  Perturbation  analysis,  and  (5)  Elasto-plastic  analysis. 

i  i  1  <  i 

;  In  order  to  develop  a  single  cpmputer  program1  using  the  selected 

'  •  ’  i 

Computational ; technique ,  it  was  nqcessary  to  formulate,  write  and 

;  '  1  i  1  I  I  .  '  i  i 

debug  programs  for  joint  perturbation  and  elasto-plastic  analysis . 

'  '  , ,  i 

The  essential  concepts  used  to  include  the  above  listed  rock 

•  ’  ,  I 

characteristics  are  discussed  in  the  subsequent  paragraphs. 

i  i  1 

i  '  .  ,  i 

I.  No  Tension  Analysis  ' 

'  ■  i  1  , 

The  basic  concept  used  in  the  combined  computer  program  for  the 

.  ;  1  !  ;  >  *  ,i 

no  tension  analysis  is  similar  to  that  developed’  by  Zienkiewicz, 

1  .  i  ;  1  ;  1  , 

et  al  (1968).  The  major  steps  iii  the  analysis  used  can  be 

, 
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summarized  as  follows: 


1.  Assign  initial  stresses  to  the  rock  mass,  and  calculate 
the  boundary  loads  required  on  the  cavity  face  to  simu- 

i  i 

late  the  creation  of  the  opening  and  other  structural 

t 

(  I 

loads  applied  to  the  system. 

I 

!  I 

i  , 

2.  Analyze  the  problem  as  a  linear  elastic  problem.  Add 

i  i 

the  induced  changes  in  stress  to  the  initial  stresses 
1  and  compute  the  principal  stresses. 

I  1 

I 

I  i 

3.  Determine  those  elements  in  which  tension  exists. 

f  ,  1 

If  the  material  is  assumed  incapable  of  sustaining 

1  i 

tension,  it  is  necessary  to  eliminate  the  tensile 

i 

stresses.  .This  is  done  by  applying  nodal  point  forces 

I 

calculated  to  eliminate  tensile  stresses. 

i  * 

I 

4.  The  elastic  analysis  is  repeated  for  the  calculated 
equivalent  nodal  point  forces;  stresses  are  determined 

/  in  the  elements.  The  check  for  tensile  stresses  is 
repeated: 

/ 

I  i 

5.  If,  at  the  end  of  step  (4),  principal  tensile  stresses 

are  still  present,  steps  (3)  and  (4)  are  repeated  until 

! 

there  is  no  appreciable  difference  in  magnitude  and  dis¬ 
tribution  of  stresses  with  further  iterations. 

I  ' 

■  I  .  '  i 

1  I 
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II .  Joint  Perturbation  Analysis 

The  stiffness  matrix  of  one-dimensional  joints  is  formulated 
according  to  the  procedure  developed  by  Goodman,  Taylor  and  Brekke 
(1968).  The  one -dimensional  joint  with  its  local  coordinate 
system  and  the  sign  convention  of  its  normal  and  shear  stresses 
are  illustrated  in  Figure  1.  The  approach  used  in  this  study 
for  the  nonlinear  analysis  of  joint  elements  is  similar  to  the 
stress  transfer  technique  proposed  by  Zienkiewicz  and  his  co¬ 
workers.  The  shear  strength  of  a  joint  element  depends  on  the 
cohesion,  c,  the  friction  angle,  <{>,  the  joint  roughness  and 
the  normal  stress  acting  across  the  joint.  Patton  (1966)  has 
shown  that  the  influence  of  joint  roughness  on  the  shear  strength 
along  a  rock  surface  can  be  taken  into  account  by  increasing  the 
friction  angle,  $,  of  the  joint  surface  by  an  amount  4>r ,  which 
is  the  average  angle  between  the  undulations  on  the  joint  sur¬ 
face  and  the  direction  of  sliding  along  the  joint.  Therefore, 
the  effective  friction  angle  d>e  of  a  rough  joint  surface  is 
given  by 

4e  a  <!>  +  <l>r  (1) 

and  the  shear  strength  of  a  joint  may  be  expressed  by 

Tf  =  c  +  oN  tan  <J>0  (2) 

where 

t£  =  shear  strength  of  the  joint 
e  =  cohesion  along  the  joint 

=  normal  stress  across  the  joint 
<(>  =  effective  friction  angle  of  the  joint  surface. 

c 
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If  the  normal  stress  across  the  joint  is  tensile,  it  is  assumed 
that  the  joint  is  incapable  of  resisting  any  shear  stress,  i.e. 
it  has  no  strength. 

The  procedure  used  to  account  for  the  nonlinear  behavior  of  the 
joint  elements  is  as  follows: 

1.  Initial  stresses  are  assigned  to  all  elements,  and 
the  boundary  loads  required  on  the  cavity  face  to 
simulate  the  creation  of  the  opening  and  other  loads 
applied  to  the  system  are  calculated. 

2.  The  problem  is  analyzed  as  a  linear  elastic  problem. 

From  the  computed  nodal  point  displacements,  the 
normal  and  tangential  displacements  across  the  joint 
are  calculated.  These  are  used  to  compute  the  changes 
in  normal  and  tangential  stresses  using  the  normal 
and  shear  stiffness  of  the  joint.  The  induced  changes 
in  normal  and  tangential  stresses  are  combined  with 
the  initial  stresses. 

3a.  Since  the  joint  is  assumed  incapable  of  sustaining  any 
shear  stress  if  a  tensile  normal  stress  is  present,  a 
check  is  made  to  see  if  a  tensile  normal  stress  exists 
across  the  joint.  If  a  tensile  normal  stress  is  present 
then  both  normal  and  tangential  stresses  are  eliminated 
and  the  equivalent  nodal  point  forces  around  the  joint 
are  calculated. 


WOODWARD-LUNDGREN  &  ASSOCIATES 

******* 


-19- 


b.  If  the  normal  stress  is  compressive,  the  shear  strength 
of  the  joint  is  calculated  by  Equation  (2).  If  the  tan¬ 
gential  stress  is  less  than  the  calculated  shear  strength, 
the  joint  remains  intact.  The  next  step  is  to  check  clo¬ 
sure  as  in  step  3(c).  If  the  tangential  stress  is  greater 
than  the  calculated  shear  strength,  the  joint  is  in  yield. 
The  excess  tangential  stress  which  is  the  difference  be¬ 
tween  the  tangential  stress  and  shear  strength  is  elimi¬ 
nated  and  replaced  by  equivalent  nodal  point  forces. 

c.  If  the  joint  is  in  compression,  and  it  undergoes  excess 
closure  beyond  the  allowable  closure  specified,  the 
nodal  point  displacements  around  the  joint  are  corrected 
to  limit  the  displacement  to  the  maximum  allowable  clo¬ 
sure.  Equivalent  normal  forces  for  the  corrections  in 
the  nodal  point  displacements  are  computed. 

4.  The  elastic  analysis  is  repeated  for  the  equivalent 
nodal  point  forces  calculated  in  Step  (3) .  The  initial 
stiffness  of  the  system  is  used  throughout  the  compu¬ 
tation.  On  the  basis  of  the  computed  stresses  and  dis¬ 
placements  the  checks  described  in  Step  (3)  are  repeated. 

5.  If  on  repeating  step  (3),  corrective  nodal  point  forces 
around  the  joints  are  still  present,  the  analysis  pro¬ 
ceeds  to  step  (4) .  This  iterative  process  is  continued 
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until  the  change  in  magnitude  and  distribution  of 
stresses  obtained  upon  iteration  is  negligible. 

Figure  4  illustrates  schematically  the  stress  transfer  process 
for  the  nonlinear  analysis  of  the  joint  systems  described  above. 

Ill .  Elasto-Plastic  Analysis 

In  developing  an  elasto-plastic  analysis.it  is  necessary  to 
define  a  yield  function  and  the  stress -strain  relations  before 
and  after  yield.  Prior  to  yield  it  is  assumed  that  linear  elas¬ 
tic  stress -strain  relations  are  applicable. 


Yield  Function 

In  the  present  study  the  yield  function  utilized  is  a  generali¬ 
zation  of  the  Mohr-Coulomb  hypothesis  suggested  by  Drucker  and 
Prager  (1952) .  The  yield  function  is  represented  by  the  following 
equation: 

f  3  a  Jj  +  /TJ  “  k  (3) 

where : 

a  and  k  *  material  constants 
*  first  stress  invariant 
J2  =  second  invariant  of  stress  deviation 


The  stress  invariants  may  be  expressed  in  terms  of  the  stress 
components  as  follows: 


(4) 
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The  yield  surface  expressed  by  Equation  (3)  for  a  >  0  is  a 
right  circular  cone  with  its  axis  equally  inclined  to  the 
coordinate  axes.  For  a  *  0  equation  (3)  reduces  to  the  Von 
Mises  yield  function;  the  yield  surface  is  a  right  circular 
cylinder.  These  two  yield  surfaces  are  illustrated  in  Figure  5. 


In  the  case  of  plane  strain,  Drucker  and  Prager  (1952)  have 
shown  that 


a 


_ tan  <fr 

(9  +  12  tar.2 


and 

k _ 3c 

(9  +  12  tan2  (jo*1 

where: 

c  =  the  cohesion  of  the  material 
$  =  angle  of  internal  friction  of  the  material 


(6) 


(7) 


Stress-Strain  Relations 

During  an  infinitesimal  increment  of  stress,  changes  in  strain 
are  assumed  to  be  composed  of  elastic  and  plastic  parts  if  the 
element  is  in  yield,  i.e. 

Ue)  =  {Ac}e  +  (Ac)  p  (8) 
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The  elastic  strain  increments  are  related  to  stress  increments 
by  the  generalized  Hooke’s  law  as 

{Aa}  ■  [D]  (Ae>e 

where  the  strain-stress  matrix  is 

Cl  -  v) 
v  (1 
0 

and  E  is  the  elastic  modulus  and  v  the 
linear  isotropic  elastic  material. 

Utilizing  the  Drucker  Prager  criteria  (equation  3)  Reyes  (1966) 
developed  elasto-plastic  Stress  strain  relations.  These  relations 
can  be  expressed  as  follows: 


{A0)  -  [“le.p. 

tic}p 

where 

“ll 

»12 

D13 

We.p.  ■ 

D21 

D22 

D23 

D31 

®32 

D33 

[D]  =  a  ”  v)  a . 


(9) 


v  o 

-  v)  0  (10) 

o  (1  -  2v) 

2  _ 

Poisson's  ratio  for  the 
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D  =  D  =  -2G  [h  +  h  (a  +  O  +  h  o  a  ] 
12  21  2  i  x  y  3xy 


D  =  D  =  -2G(h  t  +  h  a  t  ) 
i3  3i  v  i  xy  3  x  xy' 


D  =  D  =  -2G(h  r  +  h  o  t  ) 
23  32  i  xy  3  y  xy' 


where  shear  modulus  G  =  E/2  (1  +  v) 


>]/E 


a  -  9»2  §) 


(1  +  9az  £  ) 


3  v  K  k 


E  J2  (1  +  9a 


3  2  J  (1  +  9a2  £  ) 


bulk  modulus  K  = 


Since  the  elasto-plastic  stress-strain  relation  is  a  function  of 
the  elastic  constants,  the  yield  parameters  (a,  k)  and  the  stress 
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state,  it  is  necessary  to  keep  a  record  of  the  axial  stress  oz. 

For  the  plane  strain  case  the  change  in  stress, Aa,  is  related 
to  the  changes  in  stress  in  x  -  y  plane  in  the  elastic  range  by 

Aaz  =  v(Aa^  +A02)  (14) 

where  Aa1  and  A a2  are  changes  in  stress  in  two  principal  stress 
directions  in  x  -  y  plane.  In  the  plastic  range,  Drucker  and 
Prager  (1952)  have  shown  that  the  change  in  axial  stress  is 
given  by 

Aaz  =  j  (&ai  *  A 0£ )  -  j  (^ai  "  Ac^)  Sin  4>  (15) 

Computational  Procedure 

To  conduct  an  incremental  elasto-plastic  analysis  the  initial 
stress  approach  developed  by  Zienkiewicz,  et  al  (1969)  is  utilized. 
This  approach  is  consistent  with  that  described  previously  for  the 
no  tension  and  joint  perturbation  analyses.  The  load  is  applied  in 
a  series  of  increments.  The  initial  stress  process  approaches  the 
solution  of  a  nonlinear  problem  by  a  series  of  approximations. 

The  procedure  for  conducting  the  elasto-plastic  analysis  during 
a  typical  load  increment  can  be  summarized  as  follows: 

1.  For  the  applied  load  increment,  the  elastic  increments  of 

stress  (Aa’J^  and  corresponding  strain  {Ae'}^  are  determined. 
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2.  Add  {Ao*},  to  the  stresses  existing  ‘at  the  e!nd  of  the  last 

*  /  '  i  i 

I 

increment  {oQ}  to  obtain  {o’}.  The  corresponding  change>  in 
axial  stress  Ac  '  is  talcylatbd  according  to  one(  of  the  foi- 

,  Z  i  ' 

lowing  two  conditipns.  ■ 


(a)  If  f (o0),  -  k  ^  0  i.e.,  the  element  was  in  yield  at  the 


start  of  increment 


/  t 


Aoz'  ■■j.(Ao1  +  Ao2)  -  j  (Ao1  -  Ao2)  sin  <j> 


(16) 


t  i 


(b)  If  f  (oQ)  -  k  <  0  i.e.,  the  element  was  in  the  elastic 

i  i  f 

range 'at  the  start, of  increment 


[' 


Aaz'  ■  v(Aa^  +  Ao2)  1  ,  i  ( 

»  •  ’  .  :  1 
1  1  ,  l 

and  ,the  current  axial  stress  is  obtained  by 


(17) 


V  =  *  ioz' 


(18) 


Calculate  [f(o')  -  k]  using  equations  (3)  through  (7).  If  , 

}  t  1 

f(o0)  -  k  <  0  and  f(o-1)  -  k  <  0,  only  changes  in  elastic  strain 

,  < 

occur  (i.e.,1 there  is  no  yield).  No  further  computations 

are  required.  If  f(o  )  -■  k  >  6  and  f(o')  -  k1  <  0,  this 

implies  that  no  further  yielding  occurred  during  the  applied 
«  '  •'  ! 
load  increment.  The  corresponding  change  in  axial  stress  is 

i 

1  f  ,  ,  ' 

corrected  using  equation  (17),  and  oz '  recalculated  by  1 

equation  (18) .  1  ,  , 

.  ,  ;  ‘  "  ' 
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.  If,f(a')  -  k  ^  0  and  f(cQ)  -  k  ■  0  i.e.,  the  element  was 
in  yield  at  the  Start  of  increment,  calculate  {Ao}^,  using 
the  following  relation: 

i  i 

{4°>1  ■  p.Ue'}l  (19) 

I  r 

I  .  I  » 

where  [Dl  is  the  elasto-plastic  stress-strain  relation 

©  •  P  • 

expressed  by  equation  (12). 

t 

l 

i  t 

The  Excess  stresses  which  have  to  be  redistributed  are 
calculated  as  follows: 

i  1 

.  r 

{Aa'Uj  =  (Aa,>1  -  {Aa^  (20) 

f 

'  f 

i 

whereas  (Ac  ")i  is  computed  in  accordance  with  the  following 

'I  '  ' 

equation: 

i  :  ‘ 

(Aoz  )^  =  j  (Aoj"  +  Ac^")  -  j  (Aa^M  -  Aa2")  sin  4  (21) 

'  I 

f  . 

The  current  stresses,  which  are  stored,  are  computed  as 

1  ,  1  1 

follows : 

(c)  -  {o»)  -  {Aa"}1  (22) 

,  1  ■  1  .  4  i 

i 

and  current  strains 

i  '  '  i 

{ e }  -  f e_>  +  {Ae'K  (23) 

j  W  J  X  I 

<  ■ 
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4.  If  f(o’)  -  k  >  0  but  f(oQ)  -  k  <  0  i.e.,  the  element  goes 
from  the  elastic  to  the  plastic  range,  and  it  is  necessary  to 
find  the  intermediate  stress  value  at  which  yielding  commences. 
This  is  done  by  the  following  interpolation  procedure: 

Af  =  f  ({Act*  }x)  (24) 

{Aa'h  Q  =  AUa'L  (25) 

l  e .  p .  l 

where 

A  - 

(a'}  =  (a'}  -  {Aa'h  (26) 

e.p.  1  e.p. 

Use  {a'}  e.p,  to  compute  [D]  and  equation  (19)  to  calcu- 

late  {Aa}^.  Then  proceed  as  in  step  (3). 

5.  Considering  {Aa'^^  as  initial  stresses,  the  corresponding 
equilibrating  nodal  point  forces  {P}^  are  computed. 

6.  The  system  is  solved  for  the  loads  {P}^  and  (Ao*^  and 
{Ae'}2  are  determined. 

7.  Steps  2  to  6  are  repeated  until  the  nodal  forces  (P}e  reach 
sufficiently  small  values. 
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Experience  from  a  limited  number  of  solutions  performed  indicated 
that  the  rate  of  solution  convergence  would  depend  on  the  stress 
state  and  the  yield  constants  a  and  k  given.  In  general,  it  was 
found  that  if  the  Von  Mises  yield  criterion  is  used  i.e.,  a  =  0 
or  <j>  =  0,  the  solution  convergence  is  rapid,  three  to  six  itera¬ 
tions  being  necessary  in  any  increment.  However,  when  <j>  is  non¬ 
zero  the  solution  converges  very  slowly.  This  convergence  problem 
was  also  noted  by  Baker,  Sandhu  and  Shieh  (1969).  In  some  cases, 
it  may  be  found  necessary  to  use  various  numerical  schemes  to 
accelerate  convergence. 

Summary  of  Achievement  for  Phase  lb 

A  computational  procedure  using  the  initial  stress  (stress  transfer) 
approach  developed  by  Zienkiewicz  and  his  co-workers  was  formulated 
for  modelling  certain  classes  of  joints,  no  tension  characteristics 
in  the  rock  mass  and  elasto-plas tic  behavior.  This  procedure  was 
used  to  develop  a  computer  program  which  could  determine  the  mechan¬ 
ical  state  in  a  rock  mass  with  the  above-mentioned  characteristics 
in  the  vicinity  of  an  excavation.  Because  available  programs  did 
not  use  a  single  computational  technique,  it  was  necessary  to  write 
new  programs . 

ILLUSTRATIVE  PROBLEMS  (Phase  lc) 

Definition  of  Problems 

The  following  examples  were  analyzed  using  the  combined  finite 
element  computer  program  developed  in  this  study  for  the  purpose 
of  verifying  and  demonstrating  the  usage  of  the  program. 

I.  Analysis  of  a  circular  underground  opening  with  joints 
intersecting  the  opening  to  demonstrate  the  "no  tension" 
and  "joint  perturbation"  analysis. 
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II.  Elasto-plastic  analysis  of  a  thick-walled  circular  tube 
with  the  Von  Mises  yield  criterion.  A  closed  form  solu¬ 
tion  is  available  for  this  case  for  verification. 

III.  Elasto-plastic  analysis  of  a  circular  opening  with  the 
generalized  Mohr-Coulomb  yield  criterion.  The  results 
are  compared  with  those  obtained  by  Reyes  (1966)  and 
Baker,  et  al  (1969). 

IV.  Combined  no  tension,  joint  perturbation  and  elasto-plastic 
analysis  of  a  rectangular  underground  opening  to  demon¬ 
strate  the  usage  of  the  combined  computer  program. 

Results 

I.  A  Circular  Underground  Opening  with  Joints 
To  demonstrate  the  stress  transfer  technique  for  no  tension  and 
joint  perturbation  analysis,  a  circular  underground  opening 
with  two  distinct  joints,  each  inclined  22.5°  from  the  hori¬ 
zontal  axis,  was  analyzed.  Because  of  the  symmetry,  it  is 
only  necessary  to  analyze  one  quarter  of  the  full  cross-section 
as  shown  in  Figure  6.  The  finite  element  configuration  used 
in  the  analysis  is  also  shown  in  Figure  6.  The  opening  is 
assumed  to  be  in  a  uniform  stress  field  with  a  vertical  stress 
of  1000  psi  and  a  horizontal  stress  of  250  psi  corresponding 
to  Kq  =  0.25.  Material  properties  for  the  rock  and  the  joint 
characteristics  are  presented  on  Figure  6.  The  excavation  of 
the  opening  was  simulated  analytically  by  applying  boundary 
pressures  at  the  face  of  the  opening.  These  boundary  pressures 
are  equal  in  value  and  opposite  in  sign  to  the  initial  stresses 
existing  on  the  excavated  boundary.  Figures  7  and  8  show  re¬ 
spectively  the  distribution  of  major  principal  stress  for  the 
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FIG  8  -  DISTRIBUTION  OF  MAJOR  PRINCIPAL  STRESS  AROUND  A  CIRCULAR 
OPENING  (FINAL  NO  TENSION  SOLUTION,  9  CYCLES) 
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elastic  solution  indicating  the  tensile  region  developed  and 
the  final  solution  (9  cycles  of  iteration']  with  a  cracked 
region  of  lower  tensile  stress.  As  pointed  out  by  Zienkiewicz, 
those  areas  where  tensile  stresses  do  not  get  eliminated, and 
there  is  no  significant  change  in  their  magnitude  after  a 
number  of  iterations,  are  considered  as  "cracked"  or  fissured 
areas.  The  joints  remained  elastic  for  the  case  analyzed. 

II.  Elasto-Plastic  Analysis  of  a  Thick-walled  Circular  Tube 
with  the  Von  Mises  Yield  Criterion  Subject  to  Internal 
Pressure 

Prager  and  Hodge  (1951)  have  presented  a  closed  form  solution 
of  the  above  problem.  The  material  is  assumed  to  obey  the 
Von  Mises  yield  criterion,  a  special  case  of  the  Mohr-Coulomb 
criterion  expressed  by  Eq.  (3)  in  which  <f>  ■  0  or  a  ■  0  and 
k  »  c.  The  dimensions  of  the  tube,  the  material  properties 
and  the  finite  element  idealization  of  the  problem  are  shown 
in  Figure  9. 

The  circular  tube  was  subjected  to  internal  pressure  up  to 
p  =  1.39  k  in  five  unequal  increments.  The  results  of  the 
analysis  indicated  that  the  interior  surface  reached  the 
yield  surface  after  the  first  increment  of  loading  (p  >  0.75k). 
Subsequent  increase  in  pressure  caused  the  material  near  the 
interior  surface  undergoing  elasto-plastic  behavior.  For 
each  subsequent  pressure  increment,  six  cycles  of  stress 
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FIG.  9  -  FINITE  ELEMENT  MESH  FOR  AN  ELASTO-PLASTIC  ANALYSIS  OF  A 
THICK-WALLED  CIRCULAR  TUBE  (b  =  2a) 


-37- 

redistribution  were  required  for  solution  convergence,  except 
the  second  increment  for  which  only  three  cycles  were  required. 

The  results  of  the  analysis  together  with  the  closed  form  solu¬ 
tion  are  shown  in  the  form  of  stresses  and  displacements  in 
Figures  10  through  13.  Figure  10  presents  the  distribution 
of  circumferential  stress  for  the  various  pressure  increments. 

The  apex  of  these  computed  curves  indicate  the  boundary  of  the 
plastic  zone.  It  should  be  recognized  that  the  values  of  stress 
are  plotted  at  the  centroid  of  the  element.  From  Figure  10 
the  ratio  of  the  radius  of  the  plastic  zone  boundary  to  the 
internal  radius  (p/a)  is  computed  for  each  pressure  increment 

Figures  11  and  12  present  the  distribution  of  radial  and  axial 
stress  for  all  pressure  increments.  Also  shown  on  these  figures 
is  the  p/a  corresponding  to  the  distribution.  The  results  can 
be  compared  with  the  closed  form  solution  for  the  corresponding 
p/a.  The  results  are  in  excellent  agreement.  Figure  13  shows  the 
radial  displacements  on  the  interior  and  exterior  surface  as 
functions  of  the  radius  p  of  the  elastic-plastic  boundary.  Com¬ 
parison  between  the  results  obtained  from  the  finite  element 
analysis  and  those  from  the  closed  form  solution  indicates  good 
agreement. 

III.  Elasto-plastic  Analysis  of  a  Circular  Opening  with  the 
Generalized  Mohr-Coulomb  Yield  Criterion 
For  purposes  of  demonstrating  an  elasto-plastic  analysis  for 
a  material  with  the  generalized  Mohr-Coulomb  yield  criterion, 
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•  -  INCREMENT  1  p  =  0.75k 

O  -  INCREMENT  2,  ITERATION  3,  p  =  1.00k 

A  -  INCREMENT  3,  ITERATION  6,  p  =  1.18k 

□  -  INCREMENT  4,  ITERATION  6,  p  =  1.30k 

O  -  INCREMENT  5 ,  ITERATION  6, 

P  =  1.39k  -y  V* 

\  * 

p  =  RADIUS  TO  ELASTO-  Z1  A 

PLASTIC  BOUNDARY  /  /  *  X 

/v/  /  \  / 

/V  \  v 

A  /  X  . 
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*  ^ 


p/\ 


^■>s 


INCREMENT  #1- 


- -CLOSED  FORM  SOLUTION 

(PRAGER  AND  HODGE,  1951) 

- FINITE  ELEMENT  SOLUTION 


*  i 

*  •  I 


0.21— 
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NOTE:  THE  APEX  OF  THE  CURVES  REPRESENTS  THE  BOUNDARY  BETWEEN 
THE  PLASTIC  AND  ELASTIC  REGIONS.  IN  COMPARING  CURVES 
IT  IS  NECESSARY  TO  EXAMINE  CURVES  WITH  THE  SAME  p/a. 


FIG.  10  -  DISTRIBUTION  OF  CIRCUMFERENTIAL  STRESS 
FOR  A  THICK-WALLED  CIRCULAR  TUBE 
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FIG.  11  -  DISTRIBUTION  OF  RADIAL  STRESS  FOR  A 
■THICK-WALLED  CIRCULAR  TUBE 
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FIG.  12  -  DISTRIBUTION  OF  AXIAL  STRESS  FOR  A 
THICK-WALLED  CIRCULAR  TUBE 
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a  circular  opening  shown  in  Figure  14  was  analyzed.  The  defi¬ 
nition  of  the  problem  and  the  finite  element  idealization  are 
shown  in  Figure  14.  It  was  assumed  that  the  outer  boundary 
was  free.  The  excavation  of  the  opening  was  simulated  by 
applying  boundary  pressures  at  the  face  of  the  opening.  The 
boundary  pressures  which  have  vertical  and  horizontal  compo¬ 
nents  were  applied  in  7  increments.  The  percentage  of  the 
total  boundary  pressure  applied  for  each  increment  is  tabi - 
lated  in  Figure  14. 

It  was  experienced  that  the  solution  convergence  was  very 
slow  as  compared  with  the  case  for  the  thick-walled  circular 
tube  described  above.  This  convergence  problem  was  also 
noted  by  Baker,  et  al  (1969).  For  the  purpose  of  speeding 
convergence,  the  following  over-relaxation  technique  was 
adopted. 

For  every  iteration,  the  excess  stresses  (Ao"}  as  computed 
by  Eq.  (20), which  are  to  be  balanced  by  a  set  of  body  forces, 
are  multiplied  by  an  over-relaxation  factor  to  bring  the 
element  to  the  point  below  the  yield  surface.  The  over¬ 
relaxation  factor,  Ro  (p) ,  which  will  ensure  convergence 
of  the  solution, was  found  to  be  in  the  range  between  1.0 
and  1.5.  Its  value  was  selected  in  accordance  with  the 
following  procedure. 
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FIG.  14  -  FINITE  ELEMENT  MESH  FOR  AN  ELASTO- PLASTIC  ANALYSIS  OF 
A  CIRCULAR  OPENING 
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(1)  For  every  iteration  fr,  which  is  an  index  showing  where 
the  element  was  situated  outside  the  yield  surface,  is 
calculated  according  to 

fr  -  [  £  (a  * )  -  k]  /  (k  -  aJx)  (27) 

(2)  Define  the  current  state  of  stress  by  the  following 
equation 

(o)  =  {o'}  -  (Aa”}1  •  p  (28) 

p  is  to  be  determined.  Express  [f(a)  -  k]  in  terms  of  p. 

(3)  (a)  If  f  15%,  determine  p  on  the  basis  of  the  following 

equation: 

[f(o)  -  k]  -  [k  -  f(o')]  (29) 

(b)  If  15%  <  fr  <_  10%,  a  value  of  p  was  selected  such  that 

f  (o)  -  k  <  0.75  [k  -  x (a * ) J  (30) 

(c)  If  fr  <  10%,  a  value  of  p  was  selected  such  that 

f (a)  -  k  <  0.5  [k  -  f(a')]  (31) 

Using  this  technique  for  increasing  the  rate  of  convergence, 
six  to  eight  iterations  were  still  needed  for  a  loading  incre¬ 
ment.  The  results  of  the  analysis  together  with  those  pre¬ 
sented  by  Reyes  (1966)  and  Baker,  et  al  (1969)  are  shown  in 
Figures  15  and  16.  Figure  15  shows  vertical  (tangential)  and 
horizontal  (radial)  stresses  along  a  horizontal  section  through 
the  opening.  All  three  methods  of  analysis  give  the  same 
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FIG.  15  -  VERTICAL  AND  HORIZONTAL  STRESSES  ALONG  HORIZONTAL 
SECTION  OF  A  CIRCULAR  OPENING 
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Fig.  16  -  DEFO'aMATION  ALONG  CAVITY  FACE  OF  A  CIRCULAR  OPENING  AS 
COMPUTED  BY  ELASTIC  AND  ELASTIC-PLASTIC  ANALYSIS 
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distribution  of  the  horizontal  stress.  Although  all  methods 
of  analysis  indicate  that  there  is  a  reduction  in  tangential 
stress  which  results  from  local  yielding  near  the  face  of  the 
opening,  and  the  general  pattern  is  the  same,  different  dis¬ 
tributions  were  obtained.  The  difference  in  the  tangential 
stress  distribution  may  be  partially  attributed  to  the  degree 
of  convergence  obtained  by  each  method  of  analysis.  However, 
the  magnitude  of  the  difference  is  significant  in  the  vicinity 
of  the  opening.  Further  work  is  being  conducted  to  explain  the 
reasons  for  this  difference  and  if  possible  to  make  the  results 
consistent.  The  vertical  and  horizontal  components  of  the 
deflection  at  the  face  of  the  opening  are  presented  in  Figure  16. 
It  may  be  seen  that  the  elasto-plastic  solution  seems  to  give 
a  large  increase  in  both  vertical  and  horizontal  deflections 
toward  the  opening  as  compared  with  the  elastic  solution.  The 
horizontal  deflection  for  the  solution  obtained  by  Reyes  (1966) 
does  not  appear  to  be  reliable.  The  comparison  of  the  deflec¬ 
tions  indicates  differences  of  the  order  of  20  to  30%  but  the 
qualitative  pattern  of  deflection  is  consistent. 

IV.  Analysis  of  a  Rectangular  Underground  Opening 
To  demonstrate  the  usage  of  the  combined  computer  program 
developed  in  which  all  three  types  of  analysis  -  no  tension, 
joint  perturbation  and  elasto-plastic  analysis,  are  considered, 
a  hypothetical  rectangular  opening  located  at  the  depth  of 
1000  ft  was  analyzed.  The  finite  element  idealization  is 
shown  in  Figure  17,  A  horizontal  joint  is  assumed  to  exist 

WOODWARD-IUNDGREN  &  ASSOCIATES 


r 


1 


FIG.  17  -  FINITE  ELEMENT  MESH  FOR  ANALYSIS  OF  A  RECTANGULAR  OPENING 
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50  ft.  above  the  crown  of  the  opening.  Initial  stresses  were 
assumed  to  increase  with  depth  in  accordance  with  the  gravity 
stress  field.  The  vertical  stress  is  assumed  equal  to  depth 

;  ;  i  l  , 

times  unit  weight  of  the  rock  (144  pcf  assumed)  and  the  hori- 

I 

zontal  stress  equal  to  v/l-v  times  the  vertical  stress.1  The  , 

i  _  ^  i 

problem  is  defined  in  Figure  17 j  The , excavation  for  the  opening  1 

I 

was  simulated  by  applying  the  boundary  pressures  in  four  incre- 

i 

I  , 

ments  as,  shown  in  Figure  17.  1  , 

i  !  i  ■  !  1 

I  ' 

The  results  of  the  analysis  are  presented  in  Figures  18  through 

'  .  / 

24.  Figures  18  and  19  show,  respectively,  the  distribution  of 

(  I 

the  normal  and  tangential  stress  along  the  horizontal  joint  for 
both  elastic  and  combined  solution.  Some  readjustment  of  stresses 

1  i  i 

i  i 

may’  be  noted  near  the  center  line  of  the  opening  rest lting  from 

,  II  '  '  I 

the  yielding  of  the  joijit  demerits.  It  should,  however,  be 

I 

recognized  that  only  a  few  joint  elements  have  yielded  and  there 
is'  no  large  movements  along  the  joint;  since  the  joint  is  intact 

I  ,  i 

■  ■  .  I 

over  the' majority  of  its  length.  Figures  20  and  21  show,  respec- 
f  1  / 
tively,,  the  distribution  of  the  ijiajof  principal  stress  in  the  1 

i  »  ' 

vicinity  of  the  opening  obtained  from  the  elastic  and  combined 
1  ,  ,  •  11 

solution;  Figures  22  and  23  are  for  the  distribution  of  the  minor 

principal  stress  obtained  from  the  elastic  and  combined  solution. 

i  '  ; 

It,  may  be  noted,  by  comparing  with  the  elastic,  solution',  that 

/  t  i 

i  i 

some  readjustment  of  stresses  occurs  in  the  vicinity  of  the 
opening.  The  magnitude  and  region  of  tensile  stresses  have 

!  i  1 

•  i  1  : 

been  reduced  by  the  stress  transfer  process'  to  simulate  the  no 

i  ’  f  t  : 

tension  rock  characteristic.  Figure  24  indicates 1  the  development 

t  '  , 

of  the  plastic  zone  in  the  vicinity  of  the  opening.  < 
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FIG.  18  -  DISTRIBUTION  OF  NORMAL  STRESS  ALONG  HORIZONTAL 
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TENSILE  REGION 


FIG.  20  -  DISTRIBUTION  OF  MAJOR  PRINCIPAL  STRESS  AROUND  A 
RECTANGULAR  OPENING  (ELASTIC  CASE) 
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FIG.  21  -  DISTRIBUTION  OF  MAJOR  PRINCIPAL  STRESS  AROUND  A 
RECTANGULAR  OPENING  (COMBINED  FINITE  ELEMENT 
SOLUTION) 
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TENSILE  REGION 


FIG.  22  -  DISTRIBUTION  OF  MINOR  PRINCIPAL  STRESS  AROUND  A 
RECTANGULAR  OPENING  (ELASTIC  CASE) 
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FIG.  23  -  DISTRIBUTION  OF  MINOR  PRINCIPAL  STRESS  AROUND  A 
RECTANGULAR  OPENING  (COMBINED  FINITE  ELEMENT 
SOLUTION) 
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PLASTIC  ZONES 


FIG.  24  -  DEVELOPMENT  OF  PLASTIC  ZONES  AROUND  A  RECTANGULAR 
OPENING  (COMBINED  FINITE  ELEMENT  SOLUTION) 
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Summary  of  Achievement  for  Phase  lc 

The  combined  computer  program  developed  in  Phase  lb  was  used  to 
analyze  certain  illustrative  problems.  The  results  of  these 
analyses  have  been  compared  with  closed  form  solutions  and 
the  solutions  obtained  by  other  researchers.  All  comparisons 
were  satisfactory  except  for  the  case  where  the  rock  was  assumed 
to  be  an  elasto-plastic  material  with  the  yield  function  depen¬ 
dent  on  the  cohesive  and  frictional  characteristics  of  the 
material.  Further  work  is  being  conducted  in  this  area.  Using 
this  combined  computer  program,  complex  geometries  and  heteroge¬ 
neous  rock  conditions  and  properties  can  be  accounted  for. 

WORK  IN  PROGRESS  ON  PHASE  2 

The  Phase  2  work  consists  of  the  analysis  of  case  histories 
to  compare  observed  and  predicted  performance. 

Data  of  geologic  conditions  and  performance  of  two  case  histories 
of  underground  openings,  Morrow  Point  Underground  Powerplant 
and  Straight  Creek  Tunnel,  are  being  collected  and  analyzed  to 
determine  their  suitability  for  analysis  using  the  developed 
combined  computer  program  described  previously.  Possibilities 
of  analyzing  well- instrumented  laboratory  models  are  also  being 
considered. 

Our  studies  indicate  that  Morrow  Point  Underground  Powerplant 
is  suitable  for  an  analysis.  A  brief  description  of  the  pro¬ 
ject  and  pertinent  geologic  conditions  are  described  in  the 
following  paragraphs. 


WOODWARD- IUNDGREN  &  ASSOCIATES 


-58- 


The  Morrow  Point  Underground  Powerplant  was  constructed  by  the 
U.S.  Bureau  of  Reclamation  on  the  Gunnison  River  some  20  miles 
east  of  Montrose,  Colorado.  The  powerplant  chamber  was  206  ft 
long  and  57  ft  wide  with  a  height  ranging  from  65  to  134  ft  and 
is  about  400  ft  below  the  ground  surface. 

The  powerplant  is  located  entirely  within  the  Pre-Cambrian 
metamorphic  rocks.  The  composition  of  the  metamorphic  rock 
consists  predominantly  of  mica  schist  and  quartz-mica  schist 
with  some  biotite  schist  and  pegmatite.  The  bedding  strikes 
nearly  normal  to  the  powerplant  alignment  and  dips  upstream  at 
angles  varying  from  15°  to  60°,  averaging  about  35°.  There 
are  two  distinct  shear  zones  intersecting  the  powerplant  area. 
The  lower  zone  strikes  N  40°  W,  dips  32°  E,  and  has  an  average 
thickness  of  2.5  feet.  The  upper  zone  has  a  strike  of  N  20°  E, 
a  dip  of  23°  E  and  an  average  thickness  of  1.5  ft.  The  orien¬ 
tations  of  the  three  major  joint  sets  that  intersect  the  cham¬ 
bers  are:  strike  N  63°  W  and  dip  82°  SW;  strike  N  36°  E,  and 
dip  80°  W;  strike  N-S  and  dip  43°  E„ 


Five  basic  rock  types  are  present  in  the  chamber.  Their  per¬ 
centages  as  exposed  at  the  rock  surface  of  the  chamber  are 
shown  below: 


Type 

I 

Biotite  Schist 

151 

Type 

II 

Mica  Schist 

201 

Type 

III 

Quartz-Mica  Schist  or  Micaceous 
Quartzite 

40% 
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Type  IV  Quartzite 
Type  V  Pegmatite 


20% 


51 

The  rock  behavior  at  the  Morrow  Point  Underground  Powerplant 
is  significant  from  the  point  of  view  of  rock  mechanics  in 
two  respects,  Dodd,  1967;  (1)  the  excessive  rock  movements 

on  the  south-side  rock  walls,  (2)  the  rock  movements  above 
the  chamber  as  measured  by  the  multiple-position  borehole 
extensometers  (MPBX) .  It  was  hypothesized  that  the  anomalous 
rock  behavior  on  the  chamber  walls  involved  a  mass  movement 
of  a  rock  wedge,  and  the  rock  movements  above  the  chamber  as 
measured  by  the  MPBX's  were  significantly  influenced  by  the 
presence  of  shear  zones  and  joints. 

It  appears  that  using  the  combined  computer  program  developed 
which  is  capable  of  modelling  certain  types  of  joints,  shear 
zones  and  other  geologic  discontinuities  and  "no  tension" 
properties,  it  may  be  possible  to  evaluate  the  performance  of 
the  Morrow  Point  Underground  Powerplant  excavation.  However, 
it  should  be  recognized  that  the  rock  movements  at  the 
Morrow  Point  Underground  Powerplant  are  of  three-dimensional 
behavior.  Some  idealizations  are  required  for  the  analysis 
using  the  computer  program  developed  for  plane  strain  cases. 
The  finite  element  idealization  is  being  prepared  and  work 
is  proceeding  towards  analyzing  the  excavation. 
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APPENDIX  A 


A  Combined  Computer  Program  Using  Finite 
Element  Techniques  for  Elasto-plastic ,  Joint 
Perturbation,  and  No  Tension  Analysis  of 
Underground  Openings  in  Rock 
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A  COMBINED  COMPUTER  PROGRAM 
USING  FINITE  ELEMENT  TECHNIQUES  FOR  ELASTO- PLASTIC, 
JOINT  PERTURBATION,  AND  NO  TENSION  ANALYSIS 
OF  UNDERGROUND  OPENINGS  IN  ROCK 


Identification 


The  program  which  consists  of  a  main  program  and  12  subrou¬ 
tines  (STIFF,  MODIFY,  QUAD,  TRISTF,  JTSTIF,  DANSOL ,  STRESS,  REDO, 
JTSTR,  EPLAST ,  INITST,  PRINST) ,  was  developed  by  Chin-Yung  Chang 
using  a  computer  program  written  by  E.  Wilson  and  modified  by 
Goodman,  Taylor  and  Brekke  (1968). 

Purpose 

The  combined  computer  program  was  developed  on  the  basis  of 
three  concepts:  elasto-plastic ,  joint  perturbation  and  no  tension 
analysis  for  calculating  stresses  and  displacements  for  plane 
strain  conditions  in  a  rock  mass  surrounding  underground  openings. 

The  rock  mass  may  consist  of  joints,  faults,  bedding  planes  and 
other  geologic  discontinuities,  and  exhibit  elasto-plastic  and 
"no  tension"  behavior. 

Sequence  of  Operation 

(a)  The  main  program  handles  the  initial  input  and  monitors  the 
calling  of  the  subroutines  in  a  specified  order  as  bhown  in 
Fig.  Al.  If  specif ied, for  the  last  iteration  of  the  last  incre¬ 
ment  in  an  analysis,  stresses  and  excess  stresses  to  be  re¬ 
distributed  for  two-dimensional  elements,  normal  and  tangential 
stresses  for  joint  elements,  nodal  point  displacements  and  yield 
functions  for  two-dimensional  elements  are  punched  onto  cards 
to  be  ustd  for  restarting  computation. 


(b)  Subroutine  STIFF  assembles  the  general  stiffness  matrix  for 
the  entire  structure,  adds  in  concentrated  loads  at  the  nodal 
points,  adds  in  loads  due  to  boundary  pressures  and  modifies 
the  stiffness  matrix  for  the  boundary  conditions. 


(c)  In  Subroutine  QUAD  the  constitutive  equations  are  formulated. 


(d)  Subroutine  TRISTF  forms  the  stiffness  matrix  for  triangular 
subelements  and  if  specified,  element  loadings  due  to  gravity 
are  generated. 


(e)  Subroutine  JTSTIF  forms  the  stiffness  matrix  for  each  joint 
element. 
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(£)  Subroutine  MODIFY  modifies  the  stiffness  matrix  for 
the  boundary  conditions. 


(g)  Subroutine  REDO  calculates  equilibrating  nodal  point 
forces  due  to  gravity  if  specified  for  tensile 
principal  stresses  if  the  elements  are  subjected  to 
tensile  stresses  greater  than  that  allowable  and  excess 
stresses  to  be  redistributed. 

(h)  Subroutine  BANSOL  solves  the  simultaneous  equations 
representing  the  structural  stiffness  matrix  and 
the  structural  load  vector  for  nodal  point  dis¬ 
placements  . 


(i)  Subroutine  STRESS  calculates  incremental  stresses  and 
strains ,  cumulates  stresses,  and  prints  stresses  for 
two-dimensional  elements. 


(j)  Subroutine  JTSTRS  calculates  and  prints  normal  and 
tangential  stresses,  normal  and  tangential  dis¬ 
placements  (cumulative  and  incremental)  and  excess 
normal  and  tangential  stresses  to  be  redistributed 
by  comparing  stress  with  strength  for  joint  elements. 
The  equilibrating  nodal  point  forces  are  also  com¬ 
puted  from  the  excess  stresses  and  stored  for  the 
next  iteration. 


(k)  Subroutine  EPLAST  calculates  yield  functions  and 
elasto-plastic  stress-strain  relation  for  those 
two-dimensional  elements  in  yield.  The  excess 
stresses  to  be  redistributed  are  computed  as  a 
difference  between  changes  in  stress  calculated 
from  the  elastic  stress-strain  relation  and  those 
calculated  from  the  elasto-plastic  stress -strain 
relation. 

(l)  Subroutine  INITST  generates  initial  stresses  under 
free-field  conditions. 


(m)  Subroutine  PRINST  calculates  magnitudes  and  directions 
of  the  principal  stresses. 
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Output 

The  data  describing  the  finite  element  configuration, 
the  material  properties  and  pressures  applied  to  the  exca¬ 
vated  face  to  simulate  excavation  for  the  opening  are  printed 
after  being  read.  Nodal  point  displacements  (incremental 
and  cumulative),  stresses  and  yield  functions  for  two-dimen¬ 
sional  elements;  normal  and  tangential  stresses,  normal  and 
tangential  displacements  (incremental  and  cumulative)  for 
joint  elements,  are  printed  after  each  increment  or  itera¬ 
tion.  If  specified,  for  the  last  iteration  of  the  last 
increment  in  an  analysis,  stresses  and  excess  stresses  to 
be  redistributed  for  two-dimensional  elements,  normal  and 
tangential  stresses  for  joint  elements,  nodal  point  dis¬ 
placements  and  yield  functions  for  two-dimensional  elements 
are  punched  onto  cards  to  be  used  for  restarting  computation. 


Input  Data  Procedure 


1st  CARD  TYPE:  FORMAT  (8A10)  (One  Card) 


Cols  2-80  Identifying  information  to  be  printed  with 
results . 


2nd  CARD  TYPE 

:  FORMAT 

(415,  2F10.2,  215,  3F10.5) 

(One  Card) 

Cols.  1-5 

NUMNP  - 

Number  of  nodal  points 
(maximum  900) 

6-10 

NUMEL  - 

Number  of  elements 
(maximum  800) 

11-15 

NUMMAT  - 

Number  of  different  materials 
(maximum  12) 

16-20 

NUMPC  - 

Number  of  pressure  cards 

21-30 

ACELX  - 

Acceleration  in  X-direction 

WOODWARD-LUNDGREN  &  ASSOCIATES 


■  'V— n  III  ni  r  M  I  |  —  II—  IT  ’  dNUMMt 


Cols.  31-40  ACELY  - 
41-45  NP 
46-50  NRES 


51-60  FRAC 

61-70  REFPR  - 
71-80  DEPTH  - 

3rd  CARD  TYPE:  FORMAT 

Cols.  1-5  NPRSNT  - 
6-10  NREAD  - 

11-15  NPUNCH  - 

16-20  NSTSRT  - 
4th  CARD  TYPE:  FORMAT 
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Acceleration  in  Y-direction 

Number  of  approximations  (increments) 

=  -1,  Residual  stresses  generated  from 
which  residual  load  is  calculated. 

=  0,  Residual  stresses  generated,  but 

residual  load  is  zero. 

c  1,  Residual  stresses  read  as  input  from 
which  residual  load  is  generated. 

*  2,  Residual  stresses  read  as  input,  but 

residual  load  is  zero. 

Percentage  of  maximum  tensile  stress 
considered  as  cracked  zone. 

Vertical  stress  at  the  reference  point. 

Y  -  ordinate  at  the  reference  point. 

(1615)  (One  Card) 

Present  loading  increment  number. 

*0,  no  data  from  previous  computation 
will  be  read  as  input. 

=1,  data  from  last  increment  are  read  as 
input 

«0,  data  will  not  be  punched  out  at 
the  last  iteration 

=1,  data  will  be  punched  out  at  the 

last  iteration  of  the  last  increment 

*>1,  stresses  in  R-0  directions  will  be 
printed 

(15,  F10.5) 


Cols..  1-5  ITN(N)  -  Number  of  iterations  for  Nt^1  increment 

6-15  PRATIO(N)  -  Percentage  of  total  pressure  applied 

for  Nl"  increment 


Repeat  for  each  loading  increment. 


5th  CARD  TYPE:  FORMAT  (1615)  (One  Card) 


Cols.  1-5 

MJ0INT  - 

Total  number  of 
(maximum  12) 

material 

types 

6-10 

MTENS  - 

Total  number  of 
sustain  tension 

material 

types 

for  joints 
that  can 
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6th  CARD  TYPE:  FORMAT  (1615)  (Omit  this  card  if  MJOINT=0 

on  5th  card  type) 

Cols.  1-5  MJNT(I)  -  Material  type  number  for  joint  elements. 

6-10  Same 
11-15  - 

7th  CARD  TYPE;  FORMAT  (1615)  (Omit  this  card  if  MTENS=0 

on  5th  card  type) 

Cols.  1-5  MNTEN(I)  -  Material  type  number  which  can  sustain 

tension. 

6-10  Same 
11-15  - 

8th  CARD  TYPE:  FORMAT  (15,  2F10.5)  (One  Card) 

Cols.  1-5  MTYPE  -  Material  type  number. 

6-15  RO (MTYPE)-  Mass  density  of  this  material  type. 

16-20  AKO (MTYPE)  -  Ratio  of  horizontal  to  vertical 

stress . 

9th  CARD  TYPE:  FORMAT  (8F10.5) 

Cols.  1-10  E(l,  MTYPE)-  Tensile  strength  for  normal  materials 

or  normal  stiffness  for  joint  materials 

11-20  E(2,  MTYPE)-  Modulus  in  compression  for  normal  mtls. 

or  shear  stiffness  for  joint  materials 

21-30  E(3,  MTYPE)-  Poisson's  ratio  for  normal  materials 

or  cohesion  for  joint  materials 

31-40  E(4,  MTYPE)-  Modulus  in  tension  for  normal  materials  or 

angle  of  friction  for  joint  mtls . (degrees) 

41-50  E(5k  MTYPE)-  Cohesion  for  normal  materials  or 

maximum  allowable  closure  (input  as 
negative)  for  joint  materials 

51-60  E(6,  MTYPE)-  Angle  of  friction  for  normal  mtls. 

(degrees) 

Repeat  8th  and  9th  card  types  for  all  material  types. 

10th  CARD  TYPE:  FORMAT  (15,  F5.0,  4F10.0)  (One  card  for  each 

nodal  point) 

Cols.  1-5  N  -  Nodal  point  number 

6-10  CODE  (N)-  Number  which  indicates  if  displacements 

or  forces  are  to  be  specified 

=0  UR  is  the  specified  X-load  and 
UZ  is  the  specified  Y-load 
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=  1  UR  is  the  specified  X-displacement 
and  UZ  is  the  specified  Y-load 

-2  UR  is  the  specified  X-load  and 

UZ  is  the  specified  Y-displacement 

=3  UR  is  the  specified  X-displacement 
and  UZ  is  the  specified  Y-displace¬ 
ment 


Cols.  11-20 

R(N) 

X-ordinate 

21-30 

Z  (N) 

Y-ordinate 

31-40 

UR(N) 

X-load  or  displacement 

41-50 

UZ  (N) 

Y-load  or  displacement 

Nodal  points  must  be  numbered  in  sequence.  If  nodal  point  numbers 
are  omitted,  those  omitted  are  generated  automatically  at  equal 
spacings,  between  those  specified.  The  first  and  last  must  be 
specified. 


11th  CARD  TYPE:  FORMAT  (1615)  (One  card  for  each  element) 

Cols.  1-5 

M 

Element  number 

6-10 

IX(M,1)  - 

Nodal  point  I 

11-15 

IX(M,2)  - 

Nodal  point  J 

16-20 

IX(M,3) - 

Nodal  point  K 

21-25 

IX(M,4)  - 

Nodal  point  L 

26-30 

IX (M, 5)- 

Material  number 

The  nodal  point  numbers  must  be  numbered  consecutively  proceeding 
counterclockwise  around  the  elements.  The  nodal  point  numbers 
for  any  element  must  not  differ  by  more  than  39.  If  element  num¬ 
bers  are  omitted,  those  missing  will  be  generated  by  incrementing 
the  element  number  and  each  nodal  point  number  (I,  J,  K  and  L) 
by  one,  and  assigning  the  same  material  number  as  the  last  element 
specified.  The  first  and  last  elements  must  be  specified. 

Triangular  elements  are  also  permissable,  and  are  identified  by 
repeating  the  last  nodal  point  number  (i.e.,  I,  J,  K,  K) .  For 
joint  elements,  nodal  points  must  be  numbered  I,  J,  K,  L  counter¬ 
clockwise  proceeding  along  length  of  joint  from  I  to  J  and  along 
length  from  K  to  L.  Nodal  points  I  and  L  (J  and  K)  have  different 
numbers  but  identical  coordinates. 

One -dimensional  elements  are  also  permissable  and  are  identified 
by  the  node  number  sequence  (I,  J,  J,  I). 


WOODWARD- LUNDGREN  &  ASSOCIATES 


iU*' 


A-9 


12th  CARD  TYPE:  FORMAT  (1615) 


Cols.  1-5  I 

6-10  J 
11-15" 

to. 

16-  20_ 
21-25"! 


Nodal  point  number  I  along  the  boundary  IJ 
where  the  boundary  pressure  is  applied. 

Nodal  point  number  J  along  the  boundary  IJ. 

Same  as  above;  two  nodal  point  numbers  for 
each  boundary  where  the  boundary  pressure  is 
applied. 


26-30J 


As  shown  in  Figure  A2 ,  nodal  points  I  and  J  must  be  ordered 
in  counterclockwise  order  about  centroid  of  element. 


13th  CARD  TYPE:  FORMAT  (1615) 


Cols.  1-5  NPCAV  -  Number  of  nodal  points  along  the  boundary 

where  the  boundary  pressures  are  applied. 

14th  CARD  TYPE:  FORMAT  (15,  3F.10.0) 


Cols.  1-5 

NPCA(M) 

6-15 

PSCA (M,l) 

16-25 

PSCA(M, 2) 

26-35 

PSCA (M, 3) 

Nodal  point  number  where  the  boundary 
pressure  is  applied. 

X  normal  stress  at  nodal  point  NPCA(M) 
Y  normal  stress  at  nodal  point  NPCA(M) 
XY  shear  stress  at  nodal  point  NPCA(M) 


As  shown  in  Figure  A2,  a  positive  normal  stress  (ax,  a  )  is 
compressive  while  a  positive  shear  stress  (t  ) forms  a 
clockwise  moment  about  centroid  of  element.  ^ 


As  shown  in  Figure  A2,  stresses  (ax,  Oy,  Txy)  at  the  nodal 

Eoint  are  converted  to  normal  and  shear  stress  on  the 
oundary.  Then  the  nodal  point  forces  are  calculated 
from  the  normal  and  shear  stress  along  the  boundary  assum¬ 
ing  linear  stress  distribution  along  the  boundary. 

12th,  13th  and  14th  card  types  are  neglected  if  NUMPC  =  0. 

15th  CARD  TYPE:  FORMAT  (15,  3E15.5)  (One  card  for  each  element) 

Cols.  1-5  N  -  Element  number 

6-20  STRS(N,1)  -  Initial  stress  ax 

21-35  STRS(N,2)  -  Initial  stress  Oy 

36-50  STRS(N,3)  -  Initial  stress  rXy 

This  card  type  is  neglected  if  NRES  <_  0 
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